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The compilation of global Landsat data-sets and the ever-lowering costs of 
computing now make it feasible to monitor the Earth’s land cover at Landsat 
resolutions of 30 m. In this article, we describe the methods to create global 
products of forest cover and cover change at Landsat resolutions. Nevertheless, 
there are many challenges in ensuring the creation of high-quality products. And 
we propose various ways in which the challenges can be overcome. Among the 
challenges are the need for atmospheric correction, incorrect calibration 
coefficients in some of the data-sets, the different phenologies between compila- 
tions, the need for terrain correction, the lack of consistent reference data for 
training and accuracy assessment, and the need for highly automated character- 
ization and change detection. We propose and evaluate the creation and use of 
surface reflectance products, improved selection of scenes to reduce phenological 
differences, terrain illumination correction, automated training selection, and the 
use of information extraction procedures robust to errors in training data along 
with several other issues. At several stages we use Moderate Resolution Spectro- 
radiometer data and products to assist our analysis. A global working prototype 
product of forest cover and forest cover change is included. 

Keywords: Landsat; land cover; forest cover change; global mapping; global 
monitoring 


1. Introduction 

Realizing the potential of digital earth is dependent on many issues. One of the more 
important is to have internally consistent data-sets to populate it. The availability of 
global data-sets from Landsat has the potential to significantly improve the 
characterization of the Earth’s land surface. This article outlines some of the 
opportunities offered by these data-sets with reference to land cover and specifically 
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to forest cover. It also analyzes some of the challenges in using these data-sets and 
how they can be overcome. 

Land cover change is one of the most important drivers of changes in the Earth 
System. Of all land cover changes, deforestation is one of the most significant 
because of the magnitude of the resultant transformations in biophysical and 
ecological properties. Forest cover change (FCC) is highly relevant to the global 
carbon cycle, changes in the hydrological cycle, an understanding of the causes of 
changes in biodiversity and in understanding the rates and causes of land use change 
(Band 1993, Fal 1995, Houghton 1998, Pandey 2002). As a consequence, a number 
of national and international programs call for routine monitoring of global land 
cover FCCs, including the US Global Change Research Program (USGCRP 1999), 
Global Observation for Forest and Fand Cover Dynamics (GOFC-GOFD) (Skole 
et al. 1998, Townshend et al. 2004), and the Global Climate Observing System 
(GCOS 2004). 

In recent years there has been increasing emphasis on the need for products 
derived from Fandsat resolution data. Requirements for such products are specified 
in many documents, including the Earth Science Data Record (ESDR) Community 
White Paper on Fand Cover/Fand Change (Masek et al. 2006) GOFC-GOFD’s Fine 
Resolution design documents (Skole et al. 1998, Townshend et al. 2004) and the 
Integrated Global Observation Strategy’s theme Integrated Global Observation of 
the Fand theme (Townshend et al. 2010). Fandsat-class resolutions are essential for 
land cover change detection because of the fine scale of many such changes especially 
those resulting from anthropogenic factors. A substantial proportion of the 
variability of land cover change has been shown to occur at resolutions below 250 
m (Townshend and Justice 1988). 

Estimating forest cover and measuring FCC are two of the more common uses of 
Fandsat data. Fandsat-class data have primarily been used at relatively local scales 
for forest change detection. Skole and Tucker (1993), Tucker and Townshend (2000), 
Steininger et al. (2001), Zhang et al. (2005), and Huang et al. (2007) are some of the 
few studies that executed wall-to-wall change detection at national scales. The Forest 
Resource Assessment of the United Nations’ Food and Agriculture Organization 
(FAO) carried out limited Fandsat-based sampling of change detection to assist the 
estimation of global tropical forest change rates between 1990 and 2000 (FAO 2001). 
DeFries et al. (2002) calculated global tropical forest change based on Advanced 
Very High Resolution Radiometer (AVHRR) data along with regional rates of 
changes estimated from Fandsat data. Although the last two studies involved use of 
selected Fandsat imagery and products, they were not used to carry out wall-to-wall 
change mapping for the entire study area. More recently Fandsat samples combined 
with wall-to-wall data-sets have been used to provide estimates of forest loss in the 
tropics (Hansen et al. 2008) and subsequently for the globe (Hansen et al. 2010). 

Previously, continental or even global-scale analysis using Fandsat data was 
generally regarded as not feasible. This was because of the absence of well-registered 
multi-temporal data-sets, variations in sensors, the need for intensive human input 
during post-processing, variations in spectral responses of forests, the efforts needed 
to create data-sets for accuracy assessment, and the very large computational and 
storage demands in carrying out the analysis. Another major impediment was high 
data costs for global data-sets, a factor which was reduced after charging for Fandsat 
data by the United States Geological Survey (USGS) was eliminated in 2008. 
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In this article, we discuss some of the challenges and strategies for the global use 
of Landsat data. In the present article ‘global’ refers to global wall-to-wall analysis of 
the Earth’s land surface. We first discuss selection of data-sets and how some of the 
existing global data-sets published by the USGS should be refined. Then the need for 
surface reflectance products is discussed requiring atmospheric correction of the 
Landsat data-sets. Reduction of terrain effects through rectification is also essential 
as is correction for variations in solar illumination geometry. The challenges of 
assembling training data for global analysis are outlined as well as the requirement 
for data for error estimation. Procedures for global classification are discussed finally. 

Computational power was once regarded as a major barrier to large area analysis 
of Landsat-class data. This is no longer the case as is demonstrated by the fact that 
the authors use a relatively modest set of Oracle 4150 servers, with 184 processing 
cores in total, which is capable of carrying out atmospheric correction of a complete 
Landsat global data-set in approximately 4 days and to carry out change detection 
using support vector machines (SVM) can be completed in approximately a week. 
Much more effort is normally required in the painstaking tasks of data-set assembly 
and performance assessment. 


2. Selection of data-sets 

The task of carrying out global analysis of Landsat data has been made 
immeasurably easier by the assembly of global data-sets by the USGS and NASA 
from existing archives. There are now several freely available Global Land Survey 
(GLS) data-sets (Table 1) (Gutman et al. 2008). These collections were created by 
identifying optimal scenes as close to the nominal date of the collection as available. 
Lor GLS 2005 out of the nearly 500,000 Landsat images that have been acquired 
around the globe during 2004-2007, the optimal 9500 scenes have been selected 
based on several criteria including acquisition date, cloud cover, gap-fill coverage, 
sensor choice, time of year and geographic uniformity (Gutman et al. 2008). As 
Table 1 shows some of the selected scenes had a somewhat different date of collection 
compared with the nominal one. Scenes were selected to minimize cloud cover and 
also to be closest to the date of maximum greenness (though see Section 4). 

The use of GLS data-sets is facilitated by the fact that they have been 
standardized: all images are orthorectified in a Universal Transverse Mercator 
projection with a WGS84 datum, they have been resampled by cubic convolution, 
and are in a GeoTILL data format. All epochs are registered to the GLS 2000 
standard, which also forms the benchmark for the ‘standard’ LIT product from the 
US Landsat archive. Although the data are freely available from the USGS the large 
number of images in each data-set does pose logistical issues in downloading them. 


Table 1 . Global land survey collections. 


Dataset 

Acquisition year range 

Number of Images 

GLS 1975 

1972-1987 

7592 

GLS 1990 

1984-1997 

7375 

GLS 2000 

1999-2003 

8756 

GLS 2005 

2003-2008 

10,273 

GLS 2010 

2007-2011 

9000 
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Complete GLS data-sets on hard discs have been made available through the 
University of Maryland’s Global Land Cover Facility (http://glcfapp.glcf.umd. 
edu : 8 08 O/order/index . j sp ) . 

The way in which the images were selected has been progressively refined. 
Overall, the GLS 2005 data-set is superior to any other GLS in terms of data quality 
(Franks and Headley 2008). Previous GLS data-sets were single-sensor data-sets; 
Multispectral Scanner (MSS) was used for GLS 1975, Landsat 5 for GLS 1990, and 
Landsat 7 for GLS 2000. GLS 2005 had a richer selection of imagery to choose from 
by expanding the available data sources to include both Landsat 5 and 7, as well as 
images from Advanced Spaceborne Thermal Emission and Reflection Radiometer 
(ASTER) and Earth Observing Mission 1 (EO-1) (the use of the last was restricted to 
islands and reefs). Moreover, previous GLS data-sets put a high priority on cloud- 
free coverage, but at the expense of obtaining optimal leaf-on seasonality. 
Consequently, for some regions such as the dry deciduous tropics the GLS 1990 
and GLS 2000 data-sets are less useful for mapping land cover conditions because 
many images were collected during leaf-off. In GLS 2005 greater weight is put on 
selecting leaf-on images at the expense of being as cloud-free as possible. The GLS 
2010 data-set, covering imagery from 2009 to 2011 has just been released in early 
2012 . 

Global Land Survey data-sets were preceded by NASA’s global Landsat data-sets 
created by the Earthsat Corporation known as GeoCover for 1975, 1990, and 2000 
(Tucker et al. 2004), and the Mid-Decadal GLS data-set jointly assembled by NASA 
and the USGS for 2005. Although GeoCover and GLS used essentially the same 
Landsat acquisitions for each location for the 1975, 1990, and 2000 epochs, there are 
differences: one key issue is that these earlier data-sets are at 28.5 m rather than the 
30-m resolution of the GLS data-sets. Given the improvements in the processing and 
selection of images for the GLS global collections it is recommended that these are 
used in preference to the earlier ones. 


3. Refinement of global data-sets 

Comprehensive global coverage of Landsat data has become much more reliable, 
since the adoption of Landsat 7’s Long-Term Acquisition Plan (Arvidson et al. 
2001). This plan coupled with the astonishing longevity of Landsat 5 means there are 
essentially no gaps in Landsat records for GLS 2000, 2005, and 2010. Figure 1 shows 
the coverage for GLS 1975, 1990, 2000, and 2005. The initial GLS 1975 did have 
substantial gaps notably in South America, and no images were available for the 
entire eastern Siberia (Figure 1). Gaps in the former have been almost entirely filled 
using images from the Brazilian space agency INPE, which made scenes available 
from its archives. Unfortunately no substitute images have been found for northeast 
Asia. 

Examination of the GLS data-sets revealed that some scenes had significant 
deficiencies relating in particular to inappropriate phenology of the selected images. 
This is especially troublesome for change detection since differences simply due to 
phenology may result in spurious identification of change. In theory the images that 
were chosen were close to peak greenness. A test was conducted for the Korean 
peninsula (Kim et al. 2011). For 2000 and 2005 most scenes were found to have 
images from very different times of years resulting in wildly improbable rates of 
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Figure 1. GLS data coverage. The large gap in coverage for South America for GLS 1975 has 
been filled by the data provided by INPE. 


deforestation using spectral change detection methodology. A tool was constructed 
to allow automated querying of the USGS Global Visualization Viewer (GloVis) 
database (glovis.usgs.gov/;. We used Moderate Resolution Spectroradiometer 
(MODIS) data to characterize the phenology of every scene and identify scenes 
acquired on dates within 75% of peak greenness. For 9 of the 11 scenes substitute 
images were found which had phenologies which met the above criterion. Figure 2 
shows the resultant substitution for part of one scene and the greatly improved 
change results. 

Subsequently we have performed a global analysis for GLS 2000 and 2005 and 
have identified many scenes where the selected images were collected at dates quite 
distant from the peak greenness. Figure 3 shows the substitutions for GLS 2000 and 
2005. Our analysis somewhat exaggerated the severity of the problem since it includes 
scenes with a low phenological amplitude (e.g. tropical rain forests), where we do not 
need to substitute for such images. 

It is unclear why the algorithm used for scene image selection for GLS failed at 
times to identify the best scenes. It may be that for some scenes additional images 
from non-US ground receiving stations may have been selected after initial selection 
from US holdings. 


4. Production of surface reflectance images 

In order to carry out a global analysis of Landsat data it is highly advantageous that 
the pixel values in all images represent the same physical values. The DNs in Landsat 
images do not directly represent any physical values and they vary depending on the 
Landsat mission and the agency generating the image. Coefficients are almost always 
provided with calibration coefficients so the DNs can be converted to top of the 
atmosphere radiance values. The latter however do not consistently represent surface 
conditions because of variable atmospheric effects and this in turn will hinder 
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Figure 2. Spurious change detection caused by images with very different phenologies in 
deciduous woodland in South Korea (Kim et al. 2011) when using GLS data-set (top row). 
The substituted scene for 2006 clearly produces a much more realistic depiction of change than 
when the leaf off image is used (bottom row). The left and middle columns are the input 
Landsat images with bands 5, 4, 3 shown in red, green and blue. In the forest change maps 
(right column), persisting forest, persisting nonforest, forest loss and forest gain are in green, 
light pink, red and cyan, respectively. 


automated global analysis. As a consequence it is highly desirable that we estimate 
surface reflectance by carrying out atmospheric correction so that we have more 
consistent imagery for mapping change. Also it helps in the development and use of 
cross-sensor algorithms using a common radiometric basis, for example, from 
MODIS and Landsat. In addition use of reflectances will allow us better to integrate 
observations with ground-based measurements of reflectance, and to facilitate use of 
canopy reflectance models to support, for example, the creation of biophysical 
products. 

Before describing the methods used to generate the surface reflectance values 
through atmospheric correction, significant problems with the original GLS 1990 
data-set’s gain and bias coefficients need to be noted. Nearly half of these images 
were processed by different agencies using different versions of Landsat processing 
systems. As a result, different rescale gain and bias values were reported (Figure 4). 
These coefficients are not compatible with recent coefficients published by the 
USGS, that is, they yield incorrect top-of-atmosphere (TO A) reflectance values if 
used with the equations provided (Chander et al. 2004, 2009). 

To obtain proper calibration coefficients these images need to be reprocessed 
from Level 0 using the current version of the USGS Landsat processing system. The 



Downloaded by [Professor Changlin WANG] at 09:02 24 August 2012 


International Journal of Digital Earth 


1 



Figure 3. Substitutions (red quadrangles) made for GLS 2000 and 2005 because of 
phenological problems in the original GLS data-sets. 

USGS is currently attempting to acquire Level 0 data and plans to reprocess them so 
that correct values can be derived and used. 

Atmospheric correction of the Landsat images is based on the algorithm 
developed for the MODIS (Vermote et al. 2002) based on the 6S radiative transfer 
code. It was adopted for use with Landsat data and was implemented as part of the 
Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) developed 
at NASA Goddard Space Flight Center (Masek et al. 2006, Eluang et al. 2009a, b). As 
an adaptation of the MODIS Adaptive Processing System (Justice et al. 2002) for 
processing Landsat data, the LEDAPS allows rapid processing of large quantities of 
Landsat images to produce surface reflectance products from the raw radiometry. It 
has been used to produce a surface reflectance record consisting of over 2000 
Landsat images over North America (Masek et al. 2006). Song et al. (2001) showed 
how use of surface reflectance can improve the change detection using Landsat data. 
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Figure 4. Scenes with problematic gain and bias values for GLS 1990. The color coding 
represents the different gain values for band 1 . Other bands have different values, but the same 
spatial pattern. 

The 6S radiative transfer code is used to compute the transmission, intrinsic 
reflectance, and spherical albedo for relevant atmospheric constituents (Vermote 
et al. 1997), and to calculate surface reflectance by compensating for atmospheric 
scattering and absorption effects on the TOA reflectance. The atmospheric 
constituents include ozone, water vapor, and aerosols. Ozone concentration was 
derived from the Total Ozone Mapping Spectrometer (TOMS) aboard the Nimbus-7, 
Meteor-3, and Earth Probe platforms as well as from NOAA’s Television Infrared 
Observation Satellite Program’s, Operational Vertical Sounder’s, ozone data when 
TOMS data were not available. Column water vapor was taken from the reanalysis 
data of NOAA National Centers for Environmental Prediction (NCEP) (available at 
http://dss.ucar.edU/datasets/ds090.0/). Digital topography (1 km GTopo30) and 
NCEP seal level pressure data were used to adjust Rayleigh scattering to local 
conditions. Aerosol optical thickness was directly derived from the Landsat image 
using the dark, dense vegetation method of Kaufman et al. (1997). 

Explicit atmospheric correction of MSS data is not currently possible because of 
the lack of shortwave infrared bands. Until recently it was also hindered by relatively 
poor knowledge of calibration coefficients but MSS calibration for Landsats 1-3 has 
recently been reworked and standardized (Markham and Fielder in press), so it should 
be much improved. Instead, radiometric matching techniques can be used to adjust 
per-band radiometry to match a corresponding image. While this will not result in 
true surface reflectance, we believe it will be sufficient for many types of change 
detection over the long 15-year-period between the GLS 1975 and GLS 1990. Such 
techniques have been developed and used to normalize satellite data-sets in many land 
cover change studies (e.g. Hall et al. 1991, Elvidge et al. 1995, Cohen et al. 1998). An 
alternative is to explicitly correct for Rayleigh scattering and gaseous absorption of 
MSS using historical climatology data (e.g. NCEP) and just assume a constant 
aerosol loading. 



Downloaded by [Professor Changlin WANG] at 09:02 24 August 2012 


International Journal of Digital Earth 


9 


Quality Assessment (QA) of every Landsat image is highly desirable because (1) 
errors could be introduced at any of the many steps between data acquisition and 
surface reflectance generation, and (2) QA results from one image are not usually 
extensible to another because each image is processed independent of other images. 
Because the MODIS on the Terra satellite and Landsat 7 are only half an hour apart 
following the same orbit, it has been possible to create an automated Landsat- 
MODIS Consistency Checking System that automatically matches Enhance The- 
matic Mapper Plus ETM + and MODIS observations and derives a set of agreement 
metrics (Gao et al. 2006, Feng et al. 2011). Since MODIS surface reflectance 
products have been assessed comprehensively (e.g. Vermote et al. 2002, Liang et al. 
2002, Kotchenova and Vermote 2007, Vermote and Kotchenova 2008), and each of 
the six Landsat spectral bands overlaps with a MODIS band, they can be used to 
assess the quality of Landsat surface reflectance products. Our procedures involve 
the comprehensive checking of every single Landsat scene. 

Most Landsat images have been found to have close agreement with MODIS 
reflectance products (Figure 5). The discrepancies between the Landsat and MODIS 
reflectance products are generally within the uncertainty allowed by instrument 
specifications - the greater of 0.5% absolute reflectance or 5% of the retrieved 
reflectance value. Where disagreements have been found, they are all explicable. The 
most common issues are cloud movement between Landsat and MODIS overpasses, 
saturation in a Thematic mapper (TM)/ETM+ image but not MODIS, incorrect 
rescaling gains in Landsat metadata and corrupted Landsat images which needed to be 
replaced. In a few cases problematic MODIS data caused the disagreement: examples 
include the saturation of the MODIS 2.2 micron Short-wave Infrared (SWIR) band 
over very bright targets during the first few months after the Terra launch. 
Comparisons have also been made with Aeronet data and IKONOS reflectance 
data, and analysis of inter-annual stability of Landsat 5/7 reflectance time series. All 
indicate the satisfactory quality of the Landsat reflectance values (Feng et al. 2012). 


BCR 



Figure 5. Typical agreement between ETM + and MODIS reflectance values for a single 
Landsat scene. 
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The orbit of Landsat 5 images is 8 days apart from that of Terra. Hence the 
viewing geometry of the images will be different and bidirectional effects will confuse 
the comparison. Hence it is recommended that the MODIS NBAR product is used 
where corrections for view angle effects have been made in assessing the Landsat 
surface reflectance values. While MODIS NBAR data can also be used to QA the 
Landsat 7 reflectance products, MODIS daily data provide a more direct comparison 
because Terra MODIS and Landsat images acquired on the same day are expected to 
have very similar values. Assuming temporal changes due to land cover change and 
inter-annual variability are relatively small within a Landsat scene, the MODIS 
NBAR may be used to QA pre-MODIS Landsat surface reflectance products. 
Correction of Landsat images for internal bidirectional effects especially in the 
tropics is desirable but has not yet been carried out by the authors for the global 
data-sets. 

The Landsat Reflectance products are initially being made available through the 
Global Land Cover Facility (www.landcover.org). 

Identifying pixels contaminated by cloud or cloud shadow is necessary to avoid 
them falsely flagging FCC or in providing incorrect retrievals of other land surface 
attributes (Huang et al. 2010). We have relied on the implementation of LEDAPS 
(Masek et al. 2006), which currently implements two cloud masks - a version of the 
Landsat Automated Cloud Cover Assessment algorithm (Irish 2000) and a more 
aggressive mask based on MODIS spectral tests (Ackerman et al. 1998). Shadows are 
calculated from the latter using the known solar geometry and an estimate of cloud 
height based on the temperature difference between known cloudy pixels and NCEP 
surface temperature. The two methods can be combined using a voting methodology 
in order to produce a cloud and cloud shadow mask. 


5. Reducing the effects of terrain 

The GLS data-sets have undergone orthorectification using procedures similar to 
those applied in the creation of the global GeoCover data-sets (Tucker et al. 2004). 
The National Imagery and Mapping Agency (now the National Geospatial Agency) 
provided geodetic control points to the Earth Satellite Corporation where they were 
used in data processing. The goal for the horizontal Root Mean Square Error 
(RMSE) of these data-sets is better than 50 m in positional accuracy. To match GLS 
2000 as closely as possible, the GLS 2005 data-set is based on the same Digital 
Elevation Model (DEM) inputs as GLS 2000. Within the USA, the USGS National 
Elevation Data-set’s digital topography was used, while Shuttle Radar Topography 
Mission (SRTM) data are used for the rest of the globe up to 60° in latitude. In far 
northern regions othorectification uses either the Canadian Digital Elevation Data 
(CDED) or Digital Terrain Elevation Data (DTED) topographic data-sets. 
Orthorectification significantly reduces the effort of georegistering images and 
makes the task of change detection much easier and more reliable. Nevertheless, the 
nominal RMSE values of 40-50 nr are substantially greater than the pixel size and 
hence errors due to relatively small amounts of misregistration will inevitably occur 
(Townshend et al. 1992). The quality of the output scenes depends upon the accuracy 
of the geometric control and topography (Gutman et al. 2008). Some scenes may 
display geodetic errors greater than 50 m especially in mountainous areas and where 
the DEM was of poor quality. 
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Another issue relating to terrain is effects due to variations in illumination. 
Variations in the solar zenith and azimuth angles can lead to spurious differences in 
land cover being detected. Hence, especially in more rugged areas, it will be 
important to make terrain corrections. One such simple characterization is based on 
normalizing for the Illumination Condition factor (IL) which describes the incident 
solar radiation on any terrain slope facet: 

IL = cos Z. cos S + sin Z. sin S* cos(O z — <E> 5 ) 

where Z, the solar zenith angle; S, the slope angle; <Ik, the solar azimuth angle; <P S , is 
the aspect angle of the incline surface. 

The correction calculates the dependency of reflectance on IL within discrete 
ranges of vegetation greenness and then transforms the image radiometry to 
eliminate the dependency (Tan et al. 2010) 

Figure 6 shows graphically the importance of such corrections in rugged areas 
where the calendar dates of the two images are different, though note in this case that 
the difference in date is quite modest. Without making topographic corrections 
spurious changes in land cover will frequently be identified. More sophisticated 
corrections could in principle be made allowing for diffuse radiation and reflection 
from surrounding terrain (e.g. Dubayah 1992). 


6. Training data 

Land cover characterization and change detection requires an effective classification 
algorithm which is in turn dependent on adequate training. To make global mapping 
at the Landsat resolutions feasible, automated methods requiring minimal human 


1989-09-0) 2on:-1iMU Change Map vita IC 



Figure 6. Illumination correction showing original sub-scenes (top row), the corrected images 
(bottom row) and the results of change detection without and with correction (right column). 
Many of the changes using the uncorrected images are clearly spurious. In the change maps 
(right column), persisting forest, persisting nonforest, forest loss and forest gain are in green, 
light pink, red and cyan, respectively. 
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inputs are critical. This section describes a training-data automation (TDA) 
procedure for deriving training data for forest change analysis. Their use in 
classification will be described in Section 7. 

Gathering comprehensive training data for the whole of the Earth’s land surface 
is clearly very challenging logistically, since collection of adequate training data is 
one of the most time-consuming components of land cover characterization and 
change detection using remote sensing. Although there are training data associated 
with several global land cover products at coarser resolutions of 250 m and coarser, 
their value for finer resolution Landsat data is less clear, nor are they usually made 
available in user-friendly forms for other users. For some countries such as the USA 
and Brazil considerable amounts of data are available, but the different definitions of 
land cover types throughout the world make such training data challenging to use 
without considerable care in cross comparisons. 

For the task of mapping FCC globally, we have adopted a somewhat different 
approach. Note that we are concerned with cover type and we are not using the term 
‘forest’ to indicate a type of land use. Specifically we first identify training data 
automatically (Huang et al. 2008). This procedure first uses local image windows to 
identify forest training samples (Figure 7). Automation is made possible by the 
following: (1) forests are one of the darkest vegetated surfaces in satellite images 
(Colwell 1974, Goward et al. 1994, Huemmrich and Goward 1997), (2) the presence 



Figure 7. Training-data automation: an initial selection of training data is carried out by the 
use of the peak in the red response. Then a forest index is used to refine training areas for 
nonforest (Huang et al. 2008). 
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of forest within each local window is determined using to the MODIS vegetation 
continuous fields (VCF) tree cover product (Hansen et al. 2002). 

The identified forest training pixels are then used to calculate a forest index (FI) 
measuring the likelihood of a pixel being forest in reality. The most unlikely forest 
pixels are identified as nonforest training pixels. Recognizing that the training pixels 
delineated in this way are mostly pure training samples, a neighborhood algorithm is 
used to locate training pixels close to class boundaries. Specifically, less pure forest 
pixels are chosen as forest training pixels if they immediately neighbor previously 
identified pure forest pixels. Similarly, less pure nonforest pixels are determined 
according to the FI value and nearness to nonforest training pixels. This method 
creates dense sets of global ‘training data’ adapted to local conditions, which are then 
used locally as input to the classifier, with each scene having its own unique training 
data-set. 

An alternative but related approach is to use forest cover data directly from other 
sources. The MODIS VCF product at 250 m resolution can be used to identify 
similar sized training areas on Landsat scenes with an estimated percentage cover 
taken from the MODIS product. In this approach we derive training data with 
percentage cover associated with each and then estimate percentage forest cover and 
the amount of FCC using techniques such as regression trees. 

These are two possibilities for automating the creation of global training data, 
but there will doubtless be many others. Both approaches will inevitably result in 
training data-sets with errors and hence there is the need to use information- 
extraction algorithms resistant to errors as described next. 


7. Characterization of forest cover and change 

The literature is replete with many different methods of land cover characterization. 
The authors’ particular focus is in change detection and specifically in FCC 
detection. A wide range of techniques exist for cover change analysis using satellite 
data. Comprehensive reviews and comparisons of these techniques have been 
provided in several publications (e.g. Gordon 1980, Singh 1989, Lunetta and Elvidge 
1998, Coppin et al. 2004, Lu et al. 2004). Many techniques have been tested with 
varying levels of success over small areas (e.g. Jha and Unni 1994, Cohen et al. 1998, 
Lyon et al. 1998, Tokola et al. 1999, Franklin et al. 2002). However, only a few were 
used to produce FCC products for large areas, and these typically required intensive 
human inputs in order to achieve satisfactory results (e.g. Skole and Tucker 1993, 
Townshend et al. 1995, Loveland et al. 2002, Huang et al. 2007). 

As pointed out in the previous section producing globally consistent and reliable 
FCC products at affordable costs requires highly automated change detection 
algorithms requiring minimal human inputs. To meet this challenge we have 
developed a method integrating TDA and SVM. 

Previously we have used decision tree algorithms to develop land cover 
classifications at global scales (e.g. Hansen et al. 1996, DeFries et al. 1998, Hansen 
et al. 2000). However, SVM has theoretical advantages and has demonstrated 
superior performance and hence provide a good solution for global FCC analysis. 
SVM is a group of relatively new machine learning algorithms designed to find 
optimal classification solutions (Vapnik 1998), which are achieved by using the 
optimal boundary between classes as the classification boundary. As shown in Figure 
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Figure 8. A hypothetical classification problem showing the optimal boundary defined by 
support vectors (circled points) and a pair of optimal hyperplanes for an arbitrary two-class 
problem with the classes represented by solid and empty points (left). Nonlinear optimal 
boundaries defined by SVM using the RBF kernel function (right) (Song 2010). 


8, there can be many class boundaries that are sufficient to completely distinguish 
between classes for a given training data-set. SVM uses support vectors and a pair of 
optimal hyperplanes to define the optimal boundary. The optimal hyperplanes are 
defined using a vector w and are located by maximizing the distance between them. 

The value of SVM compared with decision trees and neural networks has been 
demonstrated in several previous studies (e.g. Chan et al. 2001, Huang et al. 2002; 
Pal and Mather 2005). Subsequently we applied SVM for change detection and 
showed that the accuracies were close to those produced with intensive human 
interventions (Song et al. 2005). 

A key feature of SVM is its resistance to error in training data, which is especially 
important because of our need to automate change detection given the very large size 
of the global data-sets. In experiments we progressively added noise representing error 
into training data and we have consistently found that SVM is more resistant to the 
introduction of random errors than other methods such as decision trees or neural 
networks (Figure 9) 

Since training the SVM can be automated it is possible to automate the entire 
multi-temporal SVM classification approach for FCC analysis. Using the 1990 TM 
images and 2000 ETM + images, we tested the TDA-SVM method over 19 study 
areas selected from major forest biomes across the globe (Huang et al. 2008), with 
each test area being at least a complete Landsat scene. High-resolution IKONOS 
images and independently derived reliable reference data-sets were used to calculate 
accuracy estimates for the developed FCC products. Not only the overall accuracies 
of products were well above or near 90%, the class specific user’s and producer’s 
accuracies were also quite high, with over half of those values being over 90%, 30% 
being over 80%, and the remaining being over 70% (Table 2). 

At the conclusion of the previous section we referred to an alternative procedure 
for deriving training data using the MODIS-derived VCF product. In this alternative 
procedure resultant training data are labeled with % forest cover using the MODIS 
product; this approach allows their use in estimating changes in the percentage cover 
rather than a discrete change from forest to nonforest. 
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<A> 



(B) 



Figure 9. Impact of errors in training data on the performance of different classifiers. (A) 
Variation in classification accuracy for the maximum likelihood classifier (MLC), ARTMAP 
neural net (NN), decision tree (DT), support vector machines (SVM), and kernel perceptron 
(KP). (B) Graphical depiction of the impact of 30% training data error for DT and SVM. 
Persisting forest, persisting nonforest and forest loss are in green, yellow and red, respectively 
(Song 2010). 
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Table 2. Overall, user’s and producer’s accuracies (%) of TDA-SVM derived FCC products 
evaluated using independent reference data-set (N/A: accuracy was not calculated because the 
reference data did not have information for the concerned FCC class). 


Accuracy measure Persisting forest Persisting nonforest Forest loss Forest gain 


Eastern Virginia, USA, WRS2 path 15/row 34: overall accuracy — 92 . 7 
User’s accuracy 96.7 94.9 

Producer’s accuracy 96.7 86.2 


User’s accuracy 96.7 

Producer’s accuracy 97.9 

Central Mississippi, USA, WRS2 / 
User’s accuracy 95.8 

Producer’s accuracy 88.0 


96.1 
93.3 

overall 

92.2 

95.3 


Central western Washington, USA, WRS2 path 46/row 27: overall accuracy = 94.4 
User’s accuracy 93.4 93.0 

Producer’s accuracy 99.5 99.6 

Connecticut, USA, WRS2 path 13 /row 31: overall accuracy =96.5 
User’s accuracy 98.5 95.2 

Producer’s accuracy 99.7 92.2 

Eastern Paraguay, WRS2 path 224/row 78: overall accuracy = 94.2 
User’s accuracy 76.5 95.9 

Producer’s accuracy 93.0 91.3 

Western Paraguay, WRS2 path 228 1 row 75: overall accuracy =89.6 
User’s accuracy 89.8 96.4 

Producer’s accuracy 93.1 79.7 


88.2 

74.5 

88.3 

88.4 

8 

80.6 

95.9 

96.3 

82.2 

89.4 

85.7 

72.6 

80.7 

88.9 

‘curacy = 94. 4 

98.2 

100.0 

100.0 

70.6 

88.2 

80.5 

100.0 

77.6 

83.9 

N/A 

80.4 

N/A 

85.1 

N/A 

91.8 

N/A 


8. Post-processing 

Post-processing is frequently carried out to reduce errors in the resultant products. In 
the past this has often been carried out using human interpreters who manually 
modify images (e.g. Huang et al. 2007). For a global data-set such human 
intervention is impractical if applied to the whole globe. We have implemented a 
consistency and quality checking (CQC) process, which uses other data sources to 
identify areas with major discrepancies from existing sources. This process is 
designed to highlight scenes where the FCC products’ performance may be 
unsatisfactory. 

In the CQC process, the FCC products are used to calculate percent forest cover 
(PFC) within 10-km grid cells, which are then be compared with values from the 
MODIS VCF product. A good FCC product for a particular WRS2 tile should have 
PFC values consistent with those derived using an existing reliable land cover 
product (Figure 10). Inconsistencies between the two sets of PFC values are 
indicators of uncertainties in either the FCC product or the existing land cover 
product, or both. When a FCC product is found problematic through the CQC 
process, additional training data are introduced and the images are reclassified using 
the SVM to improve the FCC product. The additional training data will be derived 
by extracting from our global training data pools we have accumulated through 
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Figure 10. Percent forest cover (PFC) within 10-km grid cells for the 2000 epoch calculated 
using a TDA-SVM derived FCC product and the standard MODIS YCF tree cover product. 
The near 1 : 1 relationship between the two sets of PCF values is an indicator of the high quality 
of the TDA-SVM derived FCC product. 


previous land cover activities (DeFries and Townshend 1994, DeFries et al. 1995, 
1998, Hansen et al. 2000, 2002) and using ultra-fine high-resolution images from 
systems such as IKONOS, QuickBird, and OrbView. 

Doubtless there are alternatives to these proposed approaches but any must 
include a very high degree of automation because of the enormous number of scenes 
to process globally. 


9. Estimation of product error 

Estimation of errors in the resultant products is essential but poses particular 
challenges, because of the variability in appearance of land cover types across the 
globe and the absence of reliable test data-sets. High-resolution images with 5 m or 
less resolution can be a valuable source of information because for some cover types 
at least they are readily recognizable by visual interpretation. Forest cover and water 
bodies fall into this category. Figure 1 1 shows a web-based tool allowing rapid 
labeling of FCC on multi-temporal Landsat images using coregistered very fine 
resolution satellite data. 

Comparing products with others that have had their accuracy assessed is an 
alternative approach as previously shown in Figure 10. We can extend this approach 
by comparing the results with a synthesis of multiple products. Song et al. (2011) 
used a regression tree model to integrate forest cover from GLC 2000, MODIS VCF, 
GLCC, GlobCover, and the standard MODIS classification to create a 5-km 
product. Comparison with validated North American products indicates that this 
synthesis is substantially more accurate than any of the individual products. The 
synthesis can now be used elsewhere in the world to assess the reliability of our forest 
products. An alternative approach would be to use crowd-sourcing derived data as in 
the Geo-Wiki Project (Fritz et al. 2009). Considerable effort will be needed to apply 
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Global Land Cover Facility 
www.landcover.org 

Forest Cover Change Labeling Tool 
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Figure 1 1 . Web-based forest change labeling tool allowing rapid labeling of forest cover and 
change using fine-resolution imagery automatically coregistered to multi-temporal Landsat 
images. 


these methods comprehensively and create statistically robust error estimates 
globally and locally. 


10. Concluding comments and future implications 

The assembly of global Landsat data collections such as those of the GLS coupled 
with rapidly declining cost of computing and improved analytical methods makes 
feasible the monitoring of the whole Earth’s land surface using Landsat data. 
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Nevertheless, there are major challenges in the creation of internally consistent 
products with known statistical errors. For example many images in the GLS 
compilations are sub-optimal for forest cover discrimination because they are leaf- 
off. We have shown that it is often possible to obtain substitute images to replace 
them. Some gaps in the original GLS data-sets especially GLS 1975 have also been 
filled. Corrections for topography are also needed in many areas. 

Use of the data-sets is made much easier by the production of Landsat surface 
reflectance images. We have demonstrated how these can be created using methods 
based on MODIS algorithms. MODIS data can also be used to evaluate the 
performance of the Landsat surface reflectance images. With some explicable and 
identified exceptions the Landsat images provide a very good representation of 
surface reflectance. The incorrect calibration coefficients provided for many scenes in 
GLS 1990 have prevented the creation of reflectance images, but it will be possible to 
replace them with correct values. 

Land cover characterization is hindered by the paucity of consistent high-quality 
training data. For forest cover the approach we adopted to overcome this issue is 
automatically to identify ‘training data’ and then apply classification methods which 
are robust to errors in the training data. We have described methods for identifying 
potentially poorly characterized images by use of independently created products. In 
this way the number of images requiring more intensive post-processing is 
considerably reduced. Estimation of errors in the resultant products is hindered by 
the absence of globally representative reliable test data-sets. Some approaches to 
reduce this problem have been proposed in this article. These have to be of very high 
quality, otherwise estimated errors may conflate errors in the product and the test 
data-sets themselves. 



Figure 12. Preliminary global forest cover change (FCC) map derived using the GLS 2000 and 
2005 datasets. The three full resolution examples (each approximately 25 km across) show 
deforestation in Lucas do Rio Verde, Brazil (lower left), logging and harvest in south central 
Sweden (lower center), FCC primarily due to burning in northeastern China (lower right). 
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Using the methods described in this paper (Figure 12), we present a working 
prototype global product of forest cover and change defined as those areas with 
greater than 30% cover. Currently extensive estimation of errors is being carried out 
and once completed a digital version will be made available through the Global Land 
Cover Facility. Although reliable estimates of errors for this product remain to be 
calculated results of comparisons with products that have had error estimation (see 
Figure 10) encourage us to believe that our approach has significant merit. We 
believe this is the first global product of any land cover type at Landsat resolution, 
specifically for forest cover and FCC to be reported in the refereed literature. 

A fundamental problem with any approach using the GLS data-sets is that each 
contains only a single image, whereas it is well established that multi-temporal images 
are of considerable value for characterizing cover types because of their different 
phenologies. It is of course possible in principle to access the archive of Landsat data 
at USGS-EDC (EROS Data Center) for free to acquire additional images, but the 
task of selecting even one additional image for each scene would be considerable. 
Nevertheless, improvements in performance will likely require at least two images for 
each scene so that simple indicators of vegetation phenology can be used. 

Can our automated approach for training identification be applied to other types 
of land cover? There is reasonable expectation that will be possible for some 
categories; for example, water has a sufficiently distinctive spectral response that 
such an approach is likely to be successful. In the creation of a 250-m global land- 
water mask, for all those areas north of the coverage of SRTM at 55°N, Landsat data 
were used successfully to identify water bodies (Carroll et al. 2009) and that change 
detection is feasible. 

Changes in data policy were crucial in allowing the creation of global Landsat 
collections. Without the policy of making Landsat images freely available the GLS 
would have been prohibitively expensive for most users. In creating global sets of 
medium resolution images there is in principle no reason to rely solely on Landsat 
images, but data policies restrict their availability. Note that for GLS 2005 some data 
from ASTER have been used. However, use of other sources is often restricted by 
charging and copyright restrictions. Government policies may often limit use of data 
in quite complex ways. Considering only Chinese data, some data are openly 
available, whereas others are restricted in use. China-Brazil Earth Resources 
Satellites (CBERS) archives are now openly available to all (as are those in Brazil 
for CBERS), but open availability of data from the more recent Huan Jing (HJ) 
satellites is openly available only to Chinese users though with some restricted 
availability for others. Beijing- 1 satellite is now operated on a commercial basis. A 
liberalization of data policies in China and elsewhere would do much to improve the 
quality of global medium resolution data-sets: for example there would be a 
considerable potential for GLS 2010 to be improved by replacing images with 
deficiencies such as poor matching phenology and high cloud cover and in selecting 
images closer the nominal 2010 date. 

One of the constraints in creating global products is dealing with all the separate 
individual scenes. The Web-Enabled Landsat Data (WELD) approach of Roy et al. 
(2010) in creating seamless mosaics for the continental US potentially has much to 
offer users if the approach was applied globally. Moreover, its selection of ‘optimal’ 
pixels from multiple images can do much to reduce the effects of clouds and other 
image artifacts. 
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In conclusion the prospect of global products derived from Landsat data at last 
allows the realization of the full potential of Landsat for global land cover 
characterization and monitoring. But successful use of existing global compilations 
requires considerable understanding of their characteristics and limitations as we 
have shown. Creating reliable global products at Landsat resolutions poses major 
challenges but the possibilities of overcoming them and the potential resultant 
opportunities have been demonstrated in this article. 
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